# models for fragmentation
library(specr)
library(lme4)
library(readr)
library(broom.mixed)
library(texreg)
library(dotwhisker)
library(dplyr)

stars <- c(0.001, 0.01, 0.05, 0.1)

# load in data
epr_yearly_group_clusters_h1 <- read_csv("data/epro_data_h1.csv")
epro_year_group_clusters <- read_csv("data/epro_data_h2.csv")
epro_year_group_clusters_geo <- read_csv("data/epro_geo_data_h3.csv")



# we have several scope conditions - not defined in data preprocessing to make them more explicit here
# 1. exclude dominant and monopoly groups (6,7) bu keeps senior and junior partners
epro_year_group_clusters_geo <- epro_year_group_clusters_geo %>%  filter(status_pwrrank <= 5)
epro_year_group_clusters <- epro_year_group_clusters %>%  filter(status_pwrrank <= 5)
epr_yearly_group_clusters_h1 <- epr_yearly_group_clusters_h1 %>%  filter(status_pwrrank <= 5)

# 2. exclude statewide groups
epro_year_group_clusters_geo <- epro_year_group_clusters_geo %>% filter(geo_typename != "Statewide")
epro_year_group_clusters <- epro_year_group_clusters %>%  filter(geo_typename != "Statewide")
epr_yearly_group_clusters_h1 <- epr_yearly_group_clusters_h1 %>%  filter(geo_typename != "Statewide")


######
## H1
#####

full_model_h1 <- lmer(n_orgs ~ resource_and_agriculture + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + (1| gwgroupid) + (1|countries_gwid),
                       data = epr_yearly_group_clusters_h1)

full_model_h1_n_res <- lmer(n_orgs ~ none_one_or_both_nats_agri_char + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + (1| gwgroupid) + (1|countries_gwid),
                      data = epr_yearly_group_clusters_h1)

screenreg(list(full_model_h1, full_model_h1_n_res), stars = stars)

#####
## H2
#####

full_model_h2 <- lmer(n_orgs ~ religious_segments_bin + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       +(1|gwgroupid) + (1|countries_gwid), data = epro_year_group_clusters)

full_model_h2_n_ed <- lmer(n_orgs ~ n_ed_religions + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                      +(1|gwgroupid) + (1|countries_gwid), data = epro_year_group_clusters)

# gives OPPOSITE effect than what we would expect!
full_model_h2_herf <- lmer(n_orgs ~ hhi_rel + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                          +(1|gwgroupid) + (1|countries_gwid), data = epro_year_group_clusters)



screenreg(list(full_model_h2,full_model_h2_n_ed,full_model_h2_herf), stars = stars)

#####
## H3
######

full_model_h3 <- lmer(n_orgs ~ n_non_intersection_group_polygons + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       +(1|gwgroupid) + (1|countries_gwid), data = epro_year_group_clusters_geo)

full_model_h3_bin <- lmer(n_orgs ~ multiple_polygons_bin + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                      +(1|gwgroupid) + (1|countries_gwid), data = epro_year_group_clusters_geo)

full_model_h3_hhi <- lmer(n_orgs ~ spatial_hhi + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                          +(1|gwgroupid) + (1|countries_gwid), data = epro_year_group_clusters_geo)

screenreg(list(full_model_h3, full_model_h3_bin, full_model_h3_hhi), stars = stars)




  








